{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "mounted-place",
   "metadata": {},
   "source": [
    "# Performing Large Numbers of Calculations with Thermo in Parallel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "confirmed-watson",
   "metadata": {},
   "source": [
    "A common request is to obtain a large number of properties from Thermo at once. Thermo is not NumPy - it cannot just automatically do all of the calculations in parallel. \n",
    "\n",
    "If you have a specific property that does not require phase equilibrium calculations to obtain, it is possible to\n",
    "use the `chemicals.numba` interface to in your own numba-accelerated code.\n",
    "https://chemicals.readthedocs.io/chemicals.numba.html\n",
    "\n",
    "For those cases where lots of flashes are needed, your best bet is to brute force it - use multiprocessing (and maybe a beefy machine) to obtain the results faster. The following code sample uses `joblib` to facilitate the calculation. Note that joblib won't show any benefits on sub-second calculations. Also note that the `threading` backend of joblib will not offer any performance improvements due to the CPython GIL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "simplified-launch",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4595970727935113"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from thermo import *\n",
    "from chemicals import *\n",
    "\n",
    "constants, properties = ChemicalConstantsPackage.from_IDs(\n",
    "    ['methane', 'ethane', 'propane', 'isobutane', 'n-butane', 'isopentane', \n",
    "     'n-pentane', 'hexane', 'heptane', 'octane', 'nonane', 'nitrogen'])\n",
    "T, P = 200, 5e6\n",
    "zs = [.8, .08, .032, .00963, .0035, .0034, .0003, .0007, .0004, .00005, .00002, .07]\n",
    "eos_kwargs = dict(Tcs=constants.Tcs, Pcs=constants.Pcs, omegas=constants.omegas)\n",
    "gas = CEOSGas(SRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "liq = CEOSLiquid(SRKMIX, eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases, T=T, P=P, zs=zs)\n",
    "# Set up a two-phase flash engine, ignoring kijs\n",
    "flasher = FlashVL(constants, properties, liquid=liq, gas=gas)\n",
    "\n",
    "# Set a composition - it could be modified in the inner loop as well\n",
    "# Do a test flash\n",
    "flasher.flash(T=T, P=P, zs=zs).gas_beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dental-liberal",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_properties(T, P):\n",
    "    # This is the function that will be called in parallel\n",
    "    # note that Python floats are faster than numpy floats\n",
    "    res = flasher.flash(T=float(T), P=float(P), zs=zs)\n",
    "    return [res.rho_mass(), res.Cp_mass(), res.gas_beta]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "opening-mainland",
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import Parallel, delayed\n",
    "pts = 30\n",
    "Ts = np.linspace(200, 400, pts)\n",
    "Ps = np.linspace(1e5, 1e7, pts)\n",
    "Ts_grid, Ps_grid = np.meshgrid(Ts, Ps)\n",
    "# processed_data = Parallel(n_jobs=16)(delayed(get_properties)(T, P) for T, P in zip(Ts_grid.flat, Ps_grid.flat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "rental-george",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# Naive loop in Python\n",
    "%timeit -r 1 -n 1 processed_data = [get_properties(T, P) for T, P in zip(Ts_grid.flat, Ps_grid.flat)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "organic-forum",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30.5 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# Use the threading feature of Joblib\n",
    "# Because the calculation is CPU-bound, the threads do not improve speed and Joblib's overhead slows down the calculation\n",
    "%timeit -r 1 -n 1 processed_data = Parallel(n_jobs=16, prefer=\"threads\")(delayed(get_properties)(T, P) for T, P in zip(Ts_grid.flat, Ps_grid.flat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "exciting-inspection",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.59 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# Use the multiprocessing feature of joblib\n",
    "# We were able to improve the speed by 5x\n",
    "%timeit -r 1 -n 1 processed_data = Parallel(n_jobs=16, batch_size=30)(delayed(get_properties)(T, P) for T, P in zip(Ts_grid.flat, Ps_grid.flat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "chubby-clock",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.98 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# For small multiprocessing jobs, the slowest job can cause a significant delay\n",
    "# For longer and larger jobs the full benefit of using all cores is shown better.\n",
    "%timeit -r 1 -n 1 processed_data = Parallel(n_jobs=8, batch_size=30)(delayed(get_properties)(T, P) for T, P in zip(Ts_grid.flat, Ps_grid.flat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "forced-entertainment",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Joblib returns the data as a flat structure, but we can re-construct it into a grid\n",
    "processed_data = Parallel(n_jobs=16, batch_size=30)(delayed(get_properties)(T, P) for T, P in zip(Ts_grid.flat, Ps_grid.flat))\n",
    "phase_fractions = np.array([[processed_data[j*pts+i][2] for j in range(pts)] for i in range(pts)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "gothic-absorption",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAG2CAYAAACH2XdzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEW0lEQVR4nO3df5yNdf7/8ecZ5qf5gTDCmClammKI2FFS22hElPp8PjY2mi2tGp9kNhu1zYg0SSuVYT4l5PNhIx+kxMaUaGkxGuljI6KRjB9fmTHDGM65vn/YOTlmMOc6x5xrznncb7frdnPe13W9r9e59to5r94/rrfNMAxDAAAAdVyQrwMAAADwBpIaAADgF0hqAACAXyCpAQAAfoGkBgAA+AWSGgAA4BdIagAAgF8gqQEAAH6BpAYAAPgFkhoAAOAXSGoAAIDWrVun/v37q0WLFrLZbFq2bNllz1m7dq1uuukmhYaGqm3btpo7d+4Vj/NSSGoAAIDKysqUlJSknJycGh2/d+9e9evXT3fccYcKCgr01FNP6dFHH9Xf/va3KxzpxdlY0BIAAJzPZrNp6dKluu+++y56zDPPPKMVK1bom2++cZb99re/1fHjx7Vq1apaiLKq+j65qoU4HA799NNPioqKks1m83U4AACLMgxDJ06cUIsWLRQUdOU6OsrLy1VRUeGVugzDqPLbFhoaqtDQUI/r3rhxo1JSUlzKUlNT9dRTT3lct1kBn9T89NNPiouL83UYAIA6Yv/+/WrVqtUVqbu8vFzXxEeq6LDdK/VFRkaqtLTUpSwrK0vjx4/3uO6ioiLFxsa6lMXGxqqkpESnTp1SeHi4x9dwV8AnNVFRUZLOPaTR0dE+jgYAYFUlJSWKi4tz/m5cCRUVFSo6bNcP+QmKjvKsNajkhEPxXfZV+X3zRiuNVQV8UlPZLBcdHU1SAwC4rNoYqhAZZVNklGfXcejK/r41b95chw4dcik7dOiQoqOjfdJKI5HUAABgOXbDIbuH03jshsM7wVxEcnKyPv74Y5ey1atXKzk5+Ype91KY0g0AgMU4ZHhlc0dpaakKCgpUUFAg6dyU7YKCAhUWFkqSxo0bp6FDhzqPHzFihL7//nv96U9/0rfffqsZM2Zo0aJFGj16tNfug7tIagAAgLZs2aLOnTurc+fOkqSMjAx17txZmZmZkqSDBw86ExxJuuaaa7RixQqtXr1aSUlJ+stf/qJZs2YpNTXVJ/FLvKdGJSUliomJUXFxMWNqAAAXVRu/F5XX+GlnK68MFG7R7seA+n1jTA0AABZjNwzZPWxz8PT8uojuJwAA4BdoqQEAwGLMDPStro5AQ1IDAIDFOGTITlLjNr/ofnrttdd0ww03KDExUU8++aQCfOwzAAABqc4nNUeOHNH06dOVn5+v7du3Kz8/X19++aWvwwIAwDRfvKfGH/hF99PZs2dVXl4uSTpz5oyaNWvm44gAADCP2U/m+LylZt26derfv79atGghm82mZcuWVTkmJydHCQkJCgsLU/fu3bVp0ybnvqZNm+rpp59W69at1aJFC6WkpKhNmza1+A0AAIAV+DypKSsrU1JSknJycqrdv3DhQmVkZCgrK0tbt25VUlKSUlNTdfjwYUnSzz//rI8++kj79u3TgQMHtGHDBq1bt642vwIAAF7l8NIWaHye1Nx999168cUXNXDgwGr3T506VcOHD1daWpoSExOVm5uriIgIzZ49W5K0Zs0atW3bVo0bN1Z4eLj69et3yTE1p0+fVklJicsGAICV2P81+8nTLdD4PKm5lIqKCuXn5yslJcVZFhQUpJSUFG3cuFGSFBcXpw0bNqi8vFx2u11r165Vu3btLlpndna2YmJinFtcXNwV/x4AALjDbnhnCzSWTmqOHj0qu92u2NhYl/LY2FgVFRVJkn7961+rb9++6ty5szp27Kg2bdpowIABF61z3LhxKi4udm779++/ot8BAADUDr+Y/TRp0iRNmjSpRseGhoYqNDT0CkcEAIB53hgTE4hjaiyd1DRp0kT16tXToUOHXMoPHTqk5s2b+ygqAACuLIdsssvmcR2BxtLdTyEhIerSpYvy8vKcZQ6HQ3l5eUpOTvZhZAAAwGp83lJTWlqq3bt3Oz/v3btXBQUFaty4sVq3bq2MjAwNGzZMXbt2Vbdu3TRt2jSVlZUpLS3Nh1EDAHDlOIxzm6d1BBqfJzVbtmzRHXfc4fyckZEhSRo2bJjmzp2rQYMG6ciRI8rMzFRRUZE6deqkVatWVRk8DACAv7B7ofvJ0/PrIpsR4Ks/lpSUKCYmRsXFxYqOjvZ1OAAAi6qN34vKa/zj/5orMsqzESKlJxzqfkNRQP2++bylBgAAuKKlxhySGgAALMZh2OQwPJz95OH5dZGlZz8BAADUFC01AABYDN1P5pDUAABgMXYFye5hZ4rdS7HUJSQ1AABYjOGFMTUGY2oAAADqpoBtqcnJyVFOTo7s9kBsoAMAWBljaswJ2Jaa9PR07dixQ5s3b/Z1KAAAuLAbQV7ZAk3gfWMAAOCXArb7CQAAq3LIJoeH7Q4OBd4qSCQ1AABYDGNqzKH7CQAA+AVaagAAsBhvDPS1G3Q/AQAAHzs3psbDBS3pfgIAAKibaKkBAMBiHF5Y+4nZTwAAwOcYU2MOSQ0AABbjUBDvqTGBMTUAAMAv0FIDAIDF2A2b7IaHL9/z8Py6iKQGAACLsXthoLCd7icAAIC6iZYaAAAsxmEEyeHh7CcHs58CR05OjnJycmS3230dCgAALuh+Midgu5/S09O1Y8cObd682dehAAAALwjYlhoAAKzKIc9nLzm8E0qdQlIDAIDFeOfle4HXGRN43xgAAPglWmoAALAY76z9FHjtFiQ1AABYjEM2OeTpmBreKAwAAHyMlhpzAu8bAwAAv0RLDQAAFuOdl+8FXrsFSQ0AABbjMGxyePqemgBcpTvw0jgAAOCXaKkBAMBiHF7ofgrEl++R1AAAYDHeWaU78JKawPvGAADAL9FSAwCAxdhlk93Dl+d5en5dRFIDAIDF0P1kTuB9YwAA4JdoqQEAwGLs8rz7yO6dUOoUkhoAACyG7idzAjapycnJUU5Ojuz2QMxlAQBWxoKW5gTeN/6X9PR07dixQ5s3b/Z1KAAAWEJOTo4SEhIUFham7t27a9OmTZc8ftq0aWrXrp3Cw8MVFxen0aNHq7y8vJairSpgkxoAAKzKkE0ODzfDzTE5CxcuVEZGhrKysrR161YlJSUpNTVVhw8frvb4BQsWaOzYscrKytI///lPvfPOO1q4cKGeffZZb9wCU0hqAACwmMruJ083d0ydOlXDhw9XWlqaEhMTlZubq4iICM2ePbva4zds2KBbbrlFgwcPVkJCgu666y49+OCDl23duZJIagAA8GMlJSUu2+nTp6scU1FRofz8fKWkpDjLgoKClJKSoo0bN1Zbb48ePZSfn+9MYr7//nt9/PHH6tu375X5IjUQsAOFAQCwKodhk8PwbEp35flxcXEu5VlZWRo/frxL2dGjR2W32xUbG+tSHhsbq2+//bba+gcPHqyjR4/q1ltvlWEYOnv2rEaMGOHT7ieSGgAALMbuhVW6K8/fv3+/oqOjneWhoaEe1Vtp7dq1eumllzRjxgx1795du3fv1qhRozRx4kQ9//zzXrmGu0hqAADwY9HR0S5JTXWaNGmievXq6dChQy7lhw4dUvPmzas95/nnn9dDDz2kRx99VJLUoUMHlZWV6bHHHtNzzz2noKDaH+HCmBoAACymsvvJ062mQkJC1KVLF+Xl5f0Sg8OhvLw8JScnV3vOyZMnqyQu9erVkyQZhmHiW3uOlhoAACzGoSA5PGx3cPf8jIwMDRs2TF27dlW3bt00bdo0lZWVKS0tTZI0dOhQtWzZUtnZ2ZKk/v37a+rUqercubOz++n5559X//79nclNbSOpAQAAGjRokI4cOaLMzEwVFRWpU6dOWrVqlXPwcGFhoUvLzJ///GfZbDb9+c9/1oEDB9S0aVP1799fkyZN8tVXkM3wVRuRRZSUlCgmJkbFxcWX7XMEAASu2vi9qLzG4+vvV2hksEd1nS49o5k9lwTU7xstNQAAWIw3p3QHEpIaAAAsxvDCKt0GC1oCAADUTbTUAABgMXbZZHdzQcrq6gg0JDUAAFiMw/B8TIwjAKcB0f0EAAD8Ai01AABYjMMLA4U9Pb8uIqkBAMBiHLLJ4eGYGE/Pr4sCL40DAAB+iZYaAAAsxm7YZPdwoLCn59dFAZvU5OTkKCcnR3a73dehAADggjE15gTeN/6X9PR07dixQ5s3b/Z1KAAAwAsCtqUGAACrcsgLaz8F4EBhkhoAACzG8MLsJ4OkBgAA+BqrdJsTsGNqAACAf6GlBgAAi2H2kzkkNQAAWAzdT+YEXhoHAAD8Ei01AABYDGs/mUNSAwCAxdD9ZA7dTwAAwC/QUgMAgMXQUmMOSQ0AABZDUmMO3U8AAMAv0FIDAIDF0FJjDkkNAAAWY8jzKdmGd0KpU0hqAACwGFpqzGFMDQAA8Au01AAAYDG01JhDUgMAgMWQ1JhD9xMAAPALtNQAAGAxtNSYQ1IDAIDFGIZNhodJiafn10V0PwEAAL9ASw0AABbjkM3jl+95en5dRFIDAIDFMKbGnIDtfsrJyVFiYqJuvvlmX4cCAAC8IGCTmvT0dO3YsUObN2/2dSgAALioHCjs6RZo6H4CAMBi6H4yh6QGAACLYUq3OSQ1AADAJ/Ly8pSXl6fDhw/L4XC47Js9e7bb9ZHUAABgMYYXup+s3lLzwgsvaMKECeratauuvvpq2Wyex0tSAwCAxRiSDMPzOqwsNzdXc+fO1UMPPeS1OgN29hMAAPCdiooK9ejRw6t1ktQAAGAxlW8U9nSzskcffVQLFizwap10PwEAYDGBMPupvLxcb731ltasWaOOHTsqODjYZf/UqVPdrpOkBgAA1Lqvv/5anTp1kiR98803LvvMDhomqQEAwGIchk02P3/53meffeb1OklqAACwGMPwwuwnq09/Os+PP/4oSWrVqpVH9TBQGAAA1DqHw6EJEyYoJiZG8fHxio+PV8OGDTVx4sQqL+KrKVpqAACwmEAYKPzcc8/pnXfe0csvv6xbbrlFkvTFF19o/PjxKi8v16RJk9yuk6QGAACLCYSk5t1339WsWbM0YMAAZ1nHjh3VsmVLPfHEEyQ1AAD4g0AYKHzs2DG1b9++Snn79u117NgxU3UypgYAANS6pKQkTZ8+vUr59OnTlZSUZKpOWmoAALCYQJj99Morr6hfv35as2aNkpOTJUkbN27U/v379fHHH5uqk5YaAAAs5lxSY/Nw8/W3uLRevXpp165dGjhwoI4fP67jx4/r/vvv186dO9WzZ09TddJSAwAAfKJFixamBgRfDEkNAAAW46+zn77++mvdeOONCgoK0tdff33JYzt27Oh2/SQ1AABYjPGvzdM6rKZTp04qKipSs2bN1KlTJ9lsNhnV9JPZbDbZ7Xa36yepAQAAtWLv3r1q2rSp89/eVqOk5v7773e74tzcXDVr1szt8wAACHT+2v0UHx/v/PcPP/ygHj16qH5911Tk7Nmz2rBhg8uxNVWj2U/Lli1TSEiIYmJiarStWLFCpaWlbgcDAAD0S/+Tp5uF3XHHHdW+ZK+4uFh33HGHqTpr3P30xhtv1LjlZfHixaaCAQAAkrzQUiMT5+fk5GjKlCkqKipSUlKS3nzzTXXr1u2ixx8/flzPPfeclixZomPHjik+Pl7Tpk1T3759Lx+eYchmqxrj//t//08NGjRwO3aphknNZ599psaNG9e40pUrV6ply5amAgIAALVv4cKFysjIUG5urrp3765p06YpNTVVO3furLZRo6KiQr1791azZs20ePFitWzZUj/88IMaNmx4yetUDmmx2Wx6+OGHFRoa6txnt9v19ddfq0ePHqa+Q42Sml69erlV6a233moqGAAA4Js3Ck+dOlXDhw9XWlqapHNjY1esWKHZs2dr7NixVY6fPXu2jh07pg0bNig4OFiSlJCQcNnrxMTE/Cs+Q1FRUQoPD3fuCwkJ0a9//WsNHz7cveD/xaPZT+Xl5aqoqHApi46O9qRKAAACnjcHCpeUlLiUh4aGurSOSOdaXfLz8zVu3DhnWVBQkFJSUrRx48Zq61++fLmSk5OVnp6uDz74QE2bNtXgwYP1zDPPqF69eheNa86cOZLOJUBjxoxRRESEqe9XHbeXSTh58qRGjhypZs2aqUGDBmrUqJHLBgAArCMuLs5lMk92dnaVY44ePSq73a7Y2FiX8tjYWBUVFVVb7/fff6/FixfLbrfr448/1vPPP6+//OUvevHFF2sU19ChQ3XgwIEq5d9995327dtXozou5HZSM2bMGH366aeaOXOmQkNDNWvWLL3wwgtq0aKF5s2bZyoIAABwHsPmnU3S/v37VVxc7NzOb43xhMPhULNmzfTWW2+pS5cuGjRokJ577jnl5ubW6PyHH35YGzZsqFL+j3/8Qw8//LCpmNzufvrwww81b9483X777UpLS1PPnj3Vtm1bxcfHa/78+RoyZIipQAAAwDneHFMTHR192aEhTZo0Ub169XTo0CGX8kOHDql58+bVnnP11VcrODjYpavp+uuvV1FRkSoqKhQSEnLJa3711Ve65ZZbqpT/+te/1siRIy957sW43VJz7NgxXXvttZLO3ajKOea33nqr1q1bZyoIX8jJyVFiYqJuvvlmX4cCAIBPhYSEqEuXLsrLy3OWORwO5eXlKTk5udpzbrnlFu3evVsOh8NZtmvXLl199dWXTWikc7OfTpw4UaW8uLjY1BIJkomk5tprr3W+2rh9+/ZatGiRpHMtOJebxmUl6enp2rFjhzZv3uzrUAAAcOWDl+9lZGTo7bff1rvvvqt//vOfevzxx1VWVuacDTV06FCXrqvHH39cx44d06hRo7Rr1y6tWLFCL730ktLT02t0vdtuu03Z2dkuCYzdbld2drbpWdRudz+lpaVp27Zt6tWrl8aOHav+/ftr+vTpOnPmjKZOnWoqCAAA8AtfLJMwaNAgHTlyRJmZmSoqKlKnTp20atUq5+DhwsJCBQX90hYSFxenv/3tbxo9erQ6duyoli1batSoUXrmmWdqdL3JkyfrtttuU7t27dSzZ09J0vr161VSUqJPP/3Urdgr2Yzqlsd0ww8//KD8/Hy1bdvW1DLhvlZSUqKYmBgVFxczHR0AcFG18XtReY3Wb2UqKCLMo7ocJ8tV+NgES/++/fTTT5o+fbq2bdum8PBwdezYUSNHjnTrhb/nq3FLjcPh0JQpU7R8+XJVVFTozjvvVFZWluLj400tOgUAAC7B4ms3eUOLFi300ksvea2+Gic1kyZN0vjx45WSkqLw8HC9/vrrOnz4sGbPnu21YAAAgP+u0l2dkydPqrCwsMrLfM30/tQ4qZk3b55mzJihP/zhD5KkNWvWqF+/fpo1a5ZLHxsAAPCQN1bZtnhLz5EjR5SWlqaVK1dWu9/MDKgaZyOFhYUuq26mpKTIZrPpp59+cvuiAAAgsD311FM6fvy4/vGPfyg8PFyrVq3Su+++q+uuu07Lly83VWeNW2rOnj2rsDDXQUvBwcE6c+aMqQsDAICLsf1r87QO6/r000/1wQcfqGvXrgoKClJ8fLx69+6t6OhoZWdnq1+/fm7XWeOkxjCMKkuEl5eXa8SIEWrQoIGzbMmSJW4HAQAAzhMA3U9lZWVq1qyZJKlRo0Y6cuSIfvWrX6lDhw7aunWrqTprnNQMGzasStnvfvc7UxcFAACBrV27dtq5c6cSEhKUlJSk//qv/1JCQoJyc3N19dVXm6qzxklN5VLhAADgCguAlppRo0bp4MGDkqSsrCz16dNH8+fPV0hIiObOnWuqTrffKAwAAK6w81bZ9qgOCzu/t6dLly764Ycf9O2336p169Zq0qSJqTprNPvp/vvvV0lJSY0rHTJkiA4fPmwqIAAA4N/OnDmjNm3a6J///KezLCIiQjfddJPphEaqYUvNBx98oCNHjtSoQsMw9OGHH2rixInOAUAAAKDmDOPc5mkdVhUcHKzy8nKv11ujpMYwDP3qV7/y+sUBAEA1AmBMTXp6uiZPnqxZs2apfn3vjIapUS2fffaZ2xW3bNnS7XMAAEBg2Lx5s/Ly8vTJJ5+oQ4cOLq+Hkcy9IqZGSU2vXr3crhgAAJgUAAOFGzZsqAceeMCrdTL7CQAAi7EZ5zZP67Ca5cuX6+6771ZwcPAVeVUMK1ECAGA1hpc2ixk4cKCOHz8uSapXr57XZ0qT1AAAgFrRtGlTffnll5LOTUKy2bzbRUb3EwAAVuOnY2pGjBihe++9VzabTTabTc2bN7/osXa73e36TSU1Z8+e1dq1a7Vnzx4NHjxYUVFR+umnnxQdHa3IyEgzVQIAgEp+OqV7/Pjx+u1vf6vdu3drwIABmjNnjho2bOi1+t1Oan744Qf16dNHhYWFOn36tHr37q2oqChNnjxZp0+fVm5urteCAwAA/qV9+/Zq3769srKy9O///u+KiIjwWt1uj6kZNWqUunbtqp9//lnh4eHO8oEDByovL89rgQEAELD8dKDw+bKysrya0EgmWmrWr1+vDRs2KCQkxKU8ISFBBw4c8FpgAAAELD/tfrrS3G6pcTgc1Q7e+fHHHxUVFeWVoAAAANzldlJz1113adq0ac7PNptNpaWlysrKUt++fb0ZGwAAgaly9pOnW4Bxu/vp1VdfVZ8+fZSYmKjy8nINHjxY3333nZo0aaK//vWvVyJGAAACir++UfhiysvLFRYW5nE9brfUxMXFadu2bXruuec0evRode7cWS+//LK++uorNWvWzOOAAACA/3M4HJo4caJatmypyMhIff/995Kk559/Xu+8846pOt1Kas6cOaM2bdrou+++05AhQ/TKK69oxowZevTRR11mQgEAAA8EwOynF198UXPnztUrr7ziMvnoxhtv1KxZs0zV6VZSExwcrPLyclMXAgAAqDRv3jy99dZbGjJkiOrVq+csT0pK0rfffmuqTre7n9LT0zV58mSdPXvW1AUBAMCl2fTLuBrTm6+/xGUcOHBAbdu2rVLucDh05swZU3W6PVB48+bNysvL0yeffKIOHTqoQYMGLvuXLFliKhAAABA4EhMTtX79esXHx7uUL168WJ07dzZVp9tJTcOGDfXAAw+YuhgAAKgBP13Q8nyZmZkaNmyYDhw4IIfDoSVLlmjnzp2aN2+ePvroI1N1up3UzJkzx9SFAABADQXAG4Xvvfdeffjhh5owYYIaNGigzMxM3XTTTfrwww/Vu3dvU3WaWqUbAADAUz179tTq1au9Vp/bSc0111wjm+3iTVqV88wBAIBJAdBSs3//ftlsNrVq1UqStGnTJi1YsECJiYl67LHHTNXpdlLz1FNPuXw+c+aMvvrqK61atUpjxowxFQQAAPhFILxRePDgwXrsscf00EMPqaioSCkpKbrxxhs1f/58FRUVKTMz0+063U5qRo0aVW15Tk6OtmzZ4nYAAAAg8HzzzTfq1q2bJGnRokXq0KGD/v73v+uTTz7RiBEjTCU1br+n5mLuvvtu/e///q+3qgMAIHAFwBuFz5w5o9DQUEnSmjVrNGDAAElS+/btdfDgQVN1ei2pWbx4sRo3buyt6gAACFwBkNTccMMNys3N1fr167V69Wr16dNHkvTTTz/pqquuMlWn291PnTt3dhkobBiGioqKdOTIEc2YMcNUEL6Qk5OjnJwc2e12X4cCAEDAmTx5sgYOHKgpU6Zo2LBhSkpKkiQtX77c2S3lLreTmvvuu8/lc1BQkJo2barbb79d7du3NxWEL6Snpys9PV0lJSWKiYnxdTgAADgFwkDh22+/XUePHlVJSYkaNWrkLH/ssccUERFhqk63k5qsrCxTFwIAADUUAG8UlqR69eq5JDSSlJCQYLo+t5OarVu3Kjg4WB06dJAkffDBB5ozZ44SExM1fvx4l+XDAQCACQHwnhrp3HjcRYsWqbCwUBUVFS77tm7d6nZ9bg8U/sMf/qBdu3ZJOveivUGDBikiIkLvv/++/vSnP7kdAAAACDxvvPGG0tLSFBsbq6+++krdunXTVVddpe+//1533323qTrdTmp27dqlTp06SZLef/999erVSwsWLNDcuXOZ0g0AgBdUjqnxdLOyGTNm6K233tKbb76pkJAQ/elPf9Lq1av15JNPqri42FSdbic1hmHI4XBIOjevvG/fvpKkuLg4HT161FQQAADgPAEwpbuwsFA9evSQJIWHh+vEiROSpIceekh//etfTdXpdlLTtWtXvfjii/rv//5vff755+rXr58kae/evYqNjTUVBAAACCzNmzfXsWPHJEmtW7fWl19+KelcPmEY5jIyt5OaadOmaevWrRo5cqSee+45tW3bVtK5wT6VGRcAAPCAN7qeLN5S85vf/EbLly+XJKWlpWn06NHq3bu3Bg0apIEDB5qq0+3ZTx07dtT27durlE+ZMkX16tUzFQQAADhPAMx+euutt5zDWdLT03XVVVdpw4YNGjBggP7whz+YqtPtpOZKLBUOAAACS1BQkIKCfukw+u1vf6vf/va3HtXpdlJz4VLhvXv31g033ODRUuEAAOA8AdBSI0nHjx/Xpk2bdPjwYWerTaWhQ4e6XZ/bSc2FS4XfeOONHi8VDgAAfhEIyyR8+OGHGjJkiEpLSxUdHe2yrqTNZjOV1Lg9UPhKLBUOAAACyx//+Ef9/ve/V2lpqY4fP66ff/7ZuVXOinKX20nNlVgqHAAABJYDBw7oySefNL14ZXXcTmomT56s//qv/9Ltt9+uBx980CtLhQMAgPMEwMv3UlNTtWXLFq/W6faYmiuxVDgAAPiFv46pqXwvjST169dPY8aM0Y4dO9ShQwcFBwe7HFs5vMUdbic10rmlEvLz87Vnzx4NHjxYUVFRCgkJIakBAAAXdd9991UpmzBhQpUym80mu93udv1uJzU//PCD+vTpo8LCQp0+fVq9e/dWVFSUJk+erNOnTys3N9ftIAAAwAUs2NLiqQunbXub22NqRo0apa5du+rnn39WeHi4s3zgwIHKy8vzanAAAASkABhTcyW4ndSsX79ef/7znxUSEuJSnpCQoAMHDngtMAAA4N/y8vJ0zz33qE2bNmrTpo3uuecerVmzxnR9bic1Doej2n6uH3/8UVFRUaYDAQAA53i6mKU3BhpfaTNmzFCfPn0UFRWlUaNGadSoUYqOjlbfvn2Vk5Njqk63x9TcddddmjZtmt566y1J5wbzlJaWKisrS3379jUVBAAAOE8ALJPw0ksv6bXXXtPIkSOdZU8++aRuueUWvfTSS0pPT3e7Trdbal599VX9/e9/V2JiosrLyzV48GBn19PkyZPdDgAAAASe48ePO1/ge7677rpLxcXFpup0u6UmLi5O27Zt08KFC7Vt2zaVlpbqkUce0ZAhQ1wGDgMAAHP89T015xswYICWLl2qMWPGuJR/8MEHuueee0zV6VZSc+bMGbVv314fffSRhgwZoiFDhpi6KAAAuAQfdT/l5ORoypQpKioqUlJSkt58880arRbw3nvv6cEHH9S9996rZcuW1ehaiYmJmjRpktauXavk5GRJ0pdffqm///3v+uMf/6g33njDeeyTTz5ZozrdSmqCg4NVXl7uzikAAKAOWLhwoTIyMpSbm6vu3btr2rRpSk1N1c6dO9WsWbOLnrdv3z49/fTT6tmzp1vXe+edd9SoUSPt2LFDO3bscJY3bNhQ77zzjvOzzWa7MkmNJKWnp2vy5MmaNWuW6tc39UJiAABwKT5oqZk6daqGDx+utLQ0SVJubq5WrFih2bNna+zYsdWeY7fbNWTIEL3wwgtav369jh8/XuPr7d27170Aa8DtrGTz5s3Ky8vTJ598og4dOqhBgwYu+5csWeK14AAACETeHFNTUlLiUh4aGqrQ0FCXsoqKCuXn52vcuHHOsqCgIKWkpGjjxo0XvcaECRPUrFkzPfLII1q/fr1nAXuB20lNw4YN9cADD1yJWAAAgOTVlpq4uDiX4qysLI0fP96l7OjRo7Lb7YqNjXUpj42N1bfffltt9V988YXeeecdFRQUeBio97id1MyZM+dKxAEAAK6A/fv3Kzo62vn5wlYaM06cOKGHHnpIb7/9tpo0aeJxfd5S46TG4XBoypQpWr58uSoqKnTnnXcqKyuLadwAAHibF1tqoqOjXZKa6jRp0kT16tXToUOHXMoPHTqk5s2bVzl+z5492rdvn/r37+8sq1yssn79+tq5c6fatGnj4RdwX41fvjdp0iQ9++yzioyMVMuWLfX666+betsfAAC4tNpeJiEkJERdunRxWZja4XAoLy/POd36fO3bt9f27dtVUFDg3AYMGKA77rhDBQUFVbq8akuNW2rmzZunGTNm6A9/+IMkac2aNerXr59mzZqloCC3X0wMAAAsJCMjQ8OGDVPXrl3VrVs3TZs2TWVlZc7ZUEOHDlXLli2VnZ2tsLAw3XjjjS7nN2zYUJKqlF/MqlWrFBkZqVtvvVXSuXfkvP3220pMTFROTo4aNWrk9neocTZSWFjosrZTSkqKbDabfvrpJ7cvCgAALsHw0uaGQYMG6dVXX1VmZqY6deqkgoICrVq1yjl4uLCwUAcPHvT8u/3LmDFjnDOztm/frj/+8Y/q27ev9u7dq4yMDFN11ril5uzZswoLC3MpCw4O1pkzZ0xdGAAAVM9XyySMHDnSZYHJ861du/aS586dO9eta+3du1eJiYmSpP/93//VPffco5deeklbt241vUB2jZMawzD08MMPu4yaLi8v14gRI1zeVcN7agAAwOWEhITo5MmTks4NaRk6dKgkqXHjxlXerVNTNU5qhg0bVqXsd7/7namLAgCAS/DR2k+16dZbb1VGRoZuueUWbdq0SQsXLpQk7dq1S61atTJVZ42TGt5PAwBALQmApGb69Ol64okntHjxYs2cOVMtW7aUJK1cuVJ9+vQxVSeLNwEAgFrXunVrffTRR1XKX3vtNdN1ktQAAGAxtn9tntZRV5SXl6uiosKl7HIvDKwOL5gBAMBqfDClu7aVlZVp5MiRatasmRo0aKBGjRq5bGaQ1AAAYDG1/UZhX/jTn/6kTz/9VDNnzlRoaKhmzZqlF154QS1atNC8efNM1Un3EwAAqHUffvih5s2bp9tvv11paWnq2bOn2rZtq/j4eM2fP19Dhgxxu05aagAAsJoA6H46duyYrr32Wknnxs8cO3ZM0rmp3uvWrTNVJ0kNAABW5McJjSRde+212rt3r6RzC2QuWrRI0rkWnMp1pNxFUgMAAGpdWlqatm3bJkkaO3ascnJyFBYWptGjR2vMmDGm6mRMDQAAFuOrtZ9q0+jRo53/TklJ0bfffqv8/Hy1bdtWHTt2NFUnSQ0AAFYTAG8UvlB8fLzi4+M9qoOkBgAA1JpTp04pLy9P99xzjyRp3LhxOn36tHN/vXr1NHHiRIWFhbldN0kNAAAW48/dT++++65WrFjhTGqmT5+uG264QeHh4ZKkb7/9Vi1atHDpnqopBgoDAGA1fjyle/78+XrsscdcyhYsWKDPPvtMn332maZMmeKcCeUukhoAAFBrdu/erQ4dOjg/h4WFKSjol3SkW7du2rFjh6m66X4CAMBi/Ln76fjx4y5jaI4cOeKy3+FwuOx3By01AABYjR93P7Vq1UrffPPNRfd//fXXatWqlam6SWoAALAaP05q+vbtq8zMTJWXl1fZd+rUKb3wwgvq16+fqbrpfgIAALXm2Wef1aJFi9SuXTuNHDlSv/rVryRJO3fu1PTp03X27Fk9++yzpuomqQEAwGL8eUxNbGysNmzYoMcff1xjx46VYZwL1GazqXfv3poxY4ZiY2NN1U1SAwCA1fj5G4WvueYarVq1SseOHdPu3bslSW3btlXjxo09qjdgk5qcnBzl5OTIbrf7OhQAAAJS48aN1a1bN6/VF7ADhdPT07Vjxw5t3rzZ16EAAODCZhhe2QJNwLbUAABgWX7e/XSlBGxLDQAA8C+01AAAYDH+PPvpSiKpAQDAauh+MoXuJwAA4BdoqQEAwGLofjKHpAYAAKuh+8kUkhoAACyGlhpzGFMDAAD8Ai01AABYDd1PppDUAABgQYHYfeQpup8AAIBfoKUGAACrMYxzm6d1BBiSGgAALIbZT+bQ/QQAAPwCLTUAAFgNs59MIakBAMBibI5zm6d1BBq6nwAAgF+gpQYAAKuh+8kUkhoAACyG2U/mkNQAAGA1vKfGFMbUAAAAv0BLDQAAFkP3kzkkNQAAWA0DhU2h+wkAAPgFWmoAALAYup/MIakBAMBqmP1kCt1PAADAL9BSAwCAxdD9ZA5JDQAAVsPsJ1PofgIAAH6BlhoAACyG7idzSGoAALAah3Fu87SOAENSAwCA1TCmxhTG1AAAAL9ASw0AABZjkxfG1HglkrqFpAYAAKvhjcKm0P0EAAD8AkkNAAAWUzml29PNXTk5OUpISFBYWJi6d++uTZs2XfTYt99+Wz179lSjRo3UqFEjpaSkXPL42kBSAwCA1Rhe2tywcOFCZWRkKCsrS1u3blVSUpJSU1N1+PDhao9fu3atHnzwQX322WfauHGj4uLidNddd+nAgQPuf18vIakBAACaOnWqhg8frrS0NCUmJio3N1cRERGaPXt2tcfPnz9fTzzxhDp16qT27dtr1qxZcjgcysvLq+XIf0FSAwCAxdgMwyubJJWUlLhsp0+frnK9iooK5efnKyUlxVkWFBSklJQUbdy4sUYxnzx5UmfOnFHjxo29cxNMIKkBAMBqHF7aJMXFxSkmJsa5ZWdnV7nc0aNHZbfbFRsb61IeGxuroqKiGoX8zDPPqEWLFi6JUW1jSjcAAH5s//79io6Odn4ODQ31+jVefvllvffee1q7dq3CwsK8Xn9NkdQAAGAx53cfeVKHJEVHR7skNdVp0qSJ6tWrp0OHDrmUHzp0SM2bN7/kua+++qpefvllrVmzRh07dvQoZk/R/QQAgNXU8uynkJAQdenSxWWQb+Wg3+Tk5Iue98orr2jixIlatWqVunbt6sYXvDJoqQEAwGp88EbhjIwMDRs2TF27dlW3bt00bdo0lZWVKS0tTZI0dOhQtWzZ0jkmZ/LkycrMzNSCBQuUkJDgHHsTGRmpyMhIz2I3iaQGAABo0KBBOnLkiDIzM1VUVKROnTpp1apVzsHDhYWFCgr6pYNn5syZqqio0L/927+51JOVlaXx48fXZuhOJDUAAFiM2TcCX1iHu0aOHKmRI0dWu2/t2rUun/ft2+f+Ba4wkhoAAKyGBS1NYaAwAADwC7TUAABgMTbHuc3TOgINSQ0AAFZD95MpdD8BAAC/QEsNAABW4+bL8y5aR4AhqQEAwGK8uUxCIKH7CQAA+AVaagAAsBoGCptCUgMAgNUYkjydkh14OQ1JDQAAVsOYGnMYUwMAAPwCLTUAAFiNIS+MqfFKJHUKSQ0AAFbDQGFT6H4CAAB+gZYaAACsxiHJ5oU6AgxJDQAAFsPsJ3PofgIAAH6BlhoAAKyGgcKmkNQAAGA1JDWm0P0EAAD8Ai01AABYDS01ppDUAABgNUzpNoWkBgAAi2FKtzmMqQEAAH6BlhoAAKyGMTWmkNQAAGA1DkOyeZiUOAIvqaH7CQAA+AVaagAAsBq6n0whqQEAwHK8kNQo8JIaup8AAIBfoKUGAACrofvJFJIaAACsxmHI4+4jZj8BAADUTbTUAABgNYbj3OZpHQGGpAYAAKthTI0pJDUAAFgNY2pMYUwNAADwC7TUAABgNXQ/mUJSAwCA1RjyQlLjlUjqFLqfAACAX6ClBgAAq6H7yRSSGgAArMbhkOThe2YcgfeeGrqfAACAX6jzSc3OnTvVqVMn5xYeHq5ly5b5OiwAAMyr7H7ydAswdb77qV27diooKJAklZaWKiEhQb179/ZtUAAAeIIxNabU+Zaa8y1fvlx33nmnGjRo4OtQAABALfN5S826des0ZcoU5efn6+DBg1q6dKnuu+8+l2NycnI0ZcoUFRUVKSkpSW+++aa6detWpa5FixZp6NChpuK4N2ao6tuCTZ0L6/rbT9t8HQIAP+E4Ya/Fi7FMghk+b6kpKytTUlKScnJyqt2/cOFCZWRkKCsrS1u3blVSUpJSU1N1+PBhl+NKSkq0YcMG9e3btzbCBgDgijEMh1e2QOPzlpq7775bd99990X3T506VcOHD1daWpokKTc3VytWrNDs2bM1duxY53EffPCB7rrrLoWFhV3yeqdPn9bp06edn0tKSjz8BgAAeJlheN7Swpgaa6moqFB+fr5SUlKcZUFBQUpJSdHGjRtdjl20aJEGDRp02Tqzs7MVExPj3OLi4rweNwAAqH2WTmqOHj0qu92u2NhYl/LY2FgVFRU5PxcXF2vTpk1KTU29bJ3jxo1TcXGxc9u/f7/X4wYAwCNM6TbF591P3hATE6NDhw7V6NjQ0FCFhoZe4YjgbQz4BRBQHA7J5uGYmAAcU2PplpomTZqoXr16VRKWQ4cOqXnz5j6KCgAAWJGlk5qQkBB16dJFeXl5zjKHw6G8vDwlJyf7MDIAAK4gup9M8Xn3U2lpqXbv3u38vHfvXhUUFKhx48Zq3bq1MjIyNGzYMHXt2lXdunXTtGnTVFZW5pwNBQCAvzEcDhkedj8xpdsHtmzZojvuuMP5OSMjQ5I0bNgwzZ07V4MGDdKRI0eUmZmpoqIiderUSatWraoyeBgAAAQ2nyc1t99+u4zLNJGNHDlSI0eOrKWIAADwMcMLbxSm+wmwJofMNaMGWXvYGABUz2FINpIad/EXHwAA+AVaagAAsBrDkEy2ULvWEVhIagAAsBjDYcjwsPvpcuNV/RFJDQAAVmM45HlLDVO6gSvu/R+/dPscuxFs6lpBNoaNAUBN5eTkaMqUKSoqKlJSUpLefPNNdevW7aLHv//++3r++ee1b98+XXfddZo8ebL69u1bixG74i8+AAAWYzgMr2zuWLhwoTIyMpSVlaWtW7cqKSlJqampOnz4cLXHb9iwQQ8++KAeeeQRffXVV7rvvvt033336ZtvvvHGLTAlYJOanJwcJSYm6uabb/Z1KAAAuDIc3tncMHXqVA0fPlxpaWlKTExUbm6uIiIiNHv27GqPf/3119WnTx+NGTNG119/vSZOnKibbrpJ06dP98YdMCVgu5/S09OVnp6u4uJiNWzYUGd1xuP3HKFmSk643897xuTrwoNtdlPnAcCFSkrP/R2qjQG43vhNOqszkqSSkhKX8tDQUIWGhrqUVVRUKD8/X+PGjXOWBQUFKSUlRRs3bqy2/o0bNzpXAaiUmpqqZcuWeRa4BwI2qal04sQJSdIX+tjHkQSO+Pa+jgAAzDtx4oRiYmKuSN0hISFq3ry5vijyzm9SZGSk4uLiXMqysrI0fvx4l7KjR4/KbrdXWYIoNjZW3377bbV1FxUVVXt8UVGR54GbFPBJTYsWLbR//35FRUXJZrM5y0tKShQXF6f9+/crOjrahxGaQ/y+Rfy+Rfy+5a/xG4ahEydOqEWLFlfs2mFhYdq7d68qKiq8Up9hGC6/bZKqtNL4k4BPaoKCgtSqVauL7o+Ojq6T/6esRPy+Rfy+Rfy+5Y/xX6kWmvOFhYUpLCzsil/nfE2aNFG9evV06NAhl/JDhw6pefPm1Z7TvHlzt46vDQE7UBgAAJwTEhKiLl26KC8vz1nmcDiUl5en5OTkas9JTk52OV6SVq9efdHja0PAt9QAAAApIyNDw4YNU9euXdWtWzdNmzZNZWVlSktLkyQNHTpULVu2VHZ2tiRp1KhR6tWrl/7yl7+oX79+eu+997Rlyxa99dZbPvsOJDUXERoaqqysrDrb90j8vkX8vkX8vkX8ddOgQYN05MgRZWZmqqioSJ06ddKqVaucg4ELCwsVFPRLB0+PHj20YMEC/fnPf9azzz6r6667TsuWLdONN97oq68gmxGIi0MAAAC/w5gaAADgF0hqAACAXyCpAQAAfoGkBgAA+IWASmqys7N18803KyoqSs2aNdN9992nnTt3uhxTXl6u9PR0XXXVVYqMjNQDDzxQ5eVChYWF6tevnyIiItSsWTONGTNGZ8+e9Xn8x44d03/+53+qXbt2Cg8PV+vWrfXkk0+quLjYpR6bzVZle++993wevyTdfvvtVWIbMWKEyzFWvf/79u2r9t7abDa9//77zuN8df9nzpypjh07Ol8olpycrJUrVzr3W/nZv1z8Vn/2Lxe/ZO1n/3LxW/3Zv9DLL78sm82mp556yllm9ecfNWQEkNTUVGPOnDnGN998YxQUFBh9+/Y1WrdubZSWljqPGTFihBEXF2fk5eUZW7ZsMX79618bPXr0cO4/e/asceONNxopKSnGV199ZXz88cdGkyZNjHHjxvk8/u3btxv333+/sXz5cmP37t1GXl6ecd111xkPPPCASz2SjDlz5hgHDx50bqdOnfJ5/IZhGL169TKGDx/uEltxcbFzv5Xv/9mzZ13iPnjwoPHCCy8YkZGRxokTJ5z1+Or+L1++3FixYoWxa9cuY+fOncazzz5rBAcHG998841hGNZ+9i8Xv9Wf/cvFbxjWfvYvF7/Vn/3zbdq0yUhISDA6duxojBo1yllu9ecfNRNQSc2FDh8+bEgyPv/8c8MwDOP48eNGcHCw8f777zuP+ec//2lIMjZu3GgYhmF8/PHHRlBQkFFUVOQ8ZubMmUZ0dLRx+vRpn8ZfnUWLFhkhISHGmTNnnGWSjKVLl9ZChJdWXfy9evVy+UNzobp2/zt16mT8/ve/dymzyv03DMNo1KiRMWvWrDr37FeqjL86Vn72K50ff1169itd6v5b8dk/ceKEcd111xmrV692ud919flHVQHV/XShyqbpxo0bS5Ly8/N15swZpaSkOI9p3769Wrdu7Vx6fePGjerQoYPLyqSpqakqKSnR//3f/9Vi9FXjv9gx0dHRql/f9T2L6enpatKkibp166bZs2fL8MHrii4W//z589WkSRPdeOONGjdunE6ePOncV5fuf35+vgoKCvTII49U2efr+2+32/Xee++prKxMycnJde7ZvzD+6lj52b9Y/HXl2b/c/bfqs5+enq5+/fq5POdS3fvbj4sL2DcKOxwOPfXUU7rlllucbz8sKipSSEiIGjZs6HLs+UupX2yp9cp9taW6+C909OhRTZw4UY899phL+YQJE/Sb3/xGERER+uSTT/TEE0+otLRUTz75ZG2ELuni8Q8ePFjx8fFq0aKFvv76az3zzDPauXOnlixZIqlu3f933nlH119/vXr06OFS7sv7v337diUnJ6u8vFyRkZFaunSpEhMTVVBQUCee/YvFfyGrPvuXir8uPPs1vf9WfPbfe+89bd26VZs3b66yry797celBWxSk56erm+++UZffPGFr0Mx5XLxl5SUqF+/fkpMTNT48eNd9j3//PPOf3fu3FllZWWaMmVKrSY1F4v//B+hDh066Oqrr9add96pPXv2qE2bNrUW3+Vc7v6fOnVKCxYscLnXlXx5/9u1a6eCggIVFxdr8eLFGjZsmD7//PMrfl1vuVj85/+wWvnZv1T8deHZr8n9t+Kzv3//fo0aNUqrV6+u9dWvUbsCsvtp5MiR+uijj/TZZ5+pVatWzvLmzZuroqJCx48fdzn+/KXUL7bUeuW+2nCx+CudOHFCffr0UVRUlJYuXarg4OBL1te9e3f9+OOPOn369JUK2cXl4r8wNknavXu3pLpx/yVp8eLFOnnypIYOHXrZ+mrz/oeEhKht27bq0qWLsrOzlZSUpNdff73OPPsXi7+S1Z/9y8V/YWyStZ79msRvxWc/Pz9fhw8f1k033aT69eurfv36+vzzz/XGG2+ofv36io2NrRPPPy4voJIawzA0cuRILV26VJ9++qmuueYal/1dunRRcHCwy1LqO3fuVGFhobPfODk5Wdu3b9fhw4edx6xevVrR0dHVNsPWZvzSuf9KveuuuxQSEqLly5fX6L9KCgoK1KhRoyu+eFtN4q8uNkm6+uqrJVn//ld65513NGDAADVt2vSy9dbW/a+Ow+HQ6dOnLf/sX0xl/JK1n/2LOT/+C1np2b+Y6uK34rN/5513avv27SooKHBuXbt21ZAhQ5z/rovPP6rhuzHKte/xxx83YmJijLVr17pMKTx58qTzmBEjRhitW7c2Pv30U2PLli1GcnKykZyc7NxfOa3vrrvuMgoKCoxVq1YZTZs2rZVpfZeLv7i42OjevbvRoUMHY/fu3S7HnD171jCMc9My3377bWP79u3Gd999Z8yYMcOIiIgwMjMzfR7/7t27jQkTJhhbtmwx9u7da3zwwQfGtddea9x2223OOqx8/yt99913hs1mM1auXFmlDl/e/7Fjxxqff/65sXfvXuPrr782xo4da9hsNuOTTz4xDMPaz/7l4rf6s3+5+K3+7F8u/kpWffarc+FsM6s//6iZgEpqJFW7zZkzx3nMqVOnjCeeeMJo1KiRERERYQwcONA4ePCgSz379u0z7r77biM8PNxo0qSJ8cc//tFl2qiv4v/ss88ueszevXsNwzCMlStXGp06dTIiIyONBg0aGElJSUZubq5ht9t9Hn9hYaFx2223GY0bNzZCQ0ONtm3bGmPGjHF5V4dhWPf+Vxo3bpwRFxdX7T315f3//e9/b8THxxshISFG06ZNjTvvvNPlB8nKz/7l4rf6s3+5+K3+7F8u/kpWffarc2FSY/XnHzVjMwwfzGcEAADwsoAaUwMAAPwXSQ0AAPALJDUAAMAvkNQAAAC/QFIDAAD8AkkNAADwCyQ1AADAL5DUAKgz1q5dK5vNJpvNpvvuu8+tcx9++GHnucuWLbsi8QHwLZIawKTKH8iLbReuEO0PEhISNG3aNF+HoZ07d2ru3LnOzw8//HCVJGfx4sUKCwvTX/7yF0nS66+/roMHD9ZilABqW31fBwDUVef/QC5cuFCZmZnauXOnsywyMtIXYbnNMAzZ7XbVr197fw4qKioUEhJi+vxmzZqpYcOGF90/a9YspaenKzc3V2lpaZKkmJgYxcTEmL4mAOujpQYwqXnz5s4tJiZGNpvNpey9997T9ddfr7CwMLVv314zZsxwnrtv3z7ZbDYtWrRIPXv2VHh4uG6++Wbt2rVLmzdvVteuXRUZGam7775bR44ccZ5X2SLxwgsvqGnTpoqOjtaIESNUUVHhPMbhcCg7O1vXXHONwsPDlZSUpMWLFzv3V3bhrFy5Ul26dFFoaKi++OIL7dmzR/fee69iY2MVGRmpm2++WWvWrHGed/vtt+uHH37Q6NGjna1RkjR+/Hh16tTJ5d5MmzZNCQkJVeKeNGmSWrRooXbt2kmS9u/fr//4j/9Qw4YN1bhxY917773at2+fR/+7vPLKK/rP//xPvffee86EBkBgoKUGuALmz5+vzMxMTZ8+XZ07d9ZXX32l4cOHq0GDBho2bJjzuKysLE2bNk2tW7fW73//ew0ePFhRUVF6/fXXFRERof/4j/9QZmamZs6c6TwnLy9PYWFhWrt2rfbt26e0tDRdddVVmjRpkiQpOztb//M//6Pc3Fxdd911WrdunX73u9+padOm6tWrl7OesWPH6tVXX9W1116rRo0aaf/+/erbt68mTZqk0NBQzZs3T/3799fOnTvVunVrLVmyRElJSXrsscc0fPhwt+9JXl6eoqOjtXr1aknSmTNnlJqaquTkZK1fv17169fXiy++qD59+ujrr7821ZLzzDPPaMaMGfroo4905513un0+gDrOxwtqAn5hzpw5RkxMjPNzmzZtjAULFrgcM3HiRCM5OdkwDMPYu3evIcmYNWuWc/9f//pXQ5KRl5fnLMvOzjbatWvn/Dxs2DCjcePGRllZmbNs5syZRmRkpGG3243y8nIjIiLC2LBhg8u1H3nkEePBBx80DOOXFa2XLVt22e91ww03GG+++abzc3x8vPHaa6+5HJOVlWUkJSW5lL322mtGfHy8S9yxsbHG6dOnnWX//d//bbRr185wOBzOstOnTxvh4eHG3/72t2rjqYz9559/dikfNmyYERISUuX+VUeSsXTp0kseA6BuoqUG8LKysjLt2bNHjzzyiEuLxtmzZ6uM6ejYsaPz37GxsZKkDh06uJQdPnzY5ZykpCRFREQ4PycnJ6u0tFT79+9XaWmpTp48qd69e7ucU1FRoc6dO7uUde3a1eVzaWmpxo8frxUrVujgwYM6e/asTp06pcLCQne+/kV16NDBpfVl27Zt2r17t6KiolyOKy8v1549e9yuv2PHjjp69KiysrLUrVu3OjOmCYD3kNQAXlZaWipJevvtt9W9e3eXffXq1XP5HBwc7Px35RiVC8scDofb116xYoVatmzpsi80NNTlc4MGDVw+P/3001q9erVeffVVtW3bVuHh4fq3f/s3l/E61QkKCpJhGC5lZ86cqXLchdcrLS1Vly5dNH/+/CrHNm3a9JLXrE7Lli21ePFi3XHHHerTp49WrlxZJWEC4N9IagAvi42NVYsWLfT9999ryJAhXq9/27ZtOnXqlMLDwyVJX375pSIjIxUXF6fGjRsrNDRUhYWFLuNnauLvf/+7Hn74YQ0cOFDSuaTjwkG7ISEhstvtLmVNmzZVUVGRDMNwJmYFBQWXvd5NN92khQsXqlmzZoqOjnYr1ouJj4/X559/7kxsVq1aRWIDBBBmPwFXwAsvvKDs7Gy98cYb2rVrl7Zv3645c+Zo6tSpHtddUVGhRx55RDt27NDHH3+srKwsjRw5UkFBQYqKitLTTz+t0aNH691339WePXu0detWvfnmm3r33XcvWe91112nJUuWqKCgQNu2bdPgwYOrtBIlJCRo3bp1OnDggI4ePSrp3KyoI0eO6JVXXtGePXuUk5OjlStXXvZ7DBkyRE2aNNG9996r9evXa+/evVq7dq2efPJJ/fjjj6bvT1xcnNauXavDhw8rNTVVJSUlpusCULeQ1ABXwKOPPqpZs2Zpzpw56tChg3r16qW5c+fqmmuu8bjuO++8U9ddd51uu+02DRo0SAMGDHB50d/EiRP1/PPPKzs7W9dff7369OmjFStWXPbaU6dOVaNGjdSjRw/1799fqampuummm1yOmTBhgvbt26c2bdo4u4iuv/56zZgxQzk5OUpKStKmTZv09NNPX/Z7REREaN26dWrdurXuv/9+XX/99XrkkUdUXl7ucctNq1attHbtWh09epTEBgggNuPCznAAlvXwww/r+PHjAfua/7Vr1+qOO+7Qzz//fMmX712KzWbT0qVL3V5mAYD10VIDoM5p1aqVHnzwQbfOGTFiBDOiAD9HSw1QhwR6S82pU6d04MABSeeWoWjevHmNzz18+LCzG+rqq6+uMhsLQN1HUgMAAPwC3U8AAMAvkNQAAAC/QFIDAAD8AkkNAADwCyQ1AADAL5DUAAAAv0BSAwAA/AJJDQAA8AskNQAAwC/8f0N7IfIfpKt+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make a plot to show the results\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "from matplotlib.colors import LogNorm\n",
    "fig, ax = plt.subplots()\n",
    "color_map = cm.viridis\n",
    "im = ax.pcolormesh(Ts_grid, Ps_grid, phase_fractions.T, cmap=color_map)\n",
    "cbar = fig.colorbar(im, ax=ax)\n",
    "cbar.set_label('Gas phase fraction')\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel('Temperature [K]')\n",
    "ax.set_ylabel('Pressure [Pa]')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}